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ABSTRACT 

Analyses and numerical procedures are presented for infrared 
radiative energy transfer in gases when other modes of energy 
transfer occur simultaneously. Two types of geometries are 
considered, a Parallel Plate Duct and a Circular Duct. Fully 
developed laminar incompressible flows of absorbing-emitting 
species in black surfaced ducts are considered under the 
conditions of uniform wall heat flux. The participating species 
considered are OH, CO, C0 2 , and H 2 0. Nongray as well as gray 
formulations are developed for both geometries. Appropriate 
limiting solutions of the governing equations are obtained and 
conduction-radiation interaction parameters are evaluated. Tien 
and Lowder's wide band model correlation has been used in 
nongray formulation. Numerical procedures are presented to 
solve the integro-differential equations for both geometries. The 
range of physical variables considered are Temperature 300 K- 
2000 K, Pressure 0.1 atm- 100.0 atm, and Spacing between 
Plates/Radius of the Tube 0.1 cm- 100 cm. An extensive 

parametric study based on nongray formulation is presented. 
Results obtained for different flow conditions indicate that the 
radiative interactions can be quite significant in fully developed 
incompressible flows. 
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Chapter 1 
INTRODUCTION 


There is renewed interest in investigating various aspects of radiative energy transfer 
in participating mediums. Radiative interactions become important in many engineering 
problems involving high temperature gases. Recent interest lies in the areas of design 
of high pressure combustion chambers and high enthalpy nozzles, entry and reentry 
phenomena, hypersonic propulsion, and defence oriented research. 

Basic formulations on radiative energy transfer in participating mediums are available 
in standard references [1-8]. The review articles presented in [9-24] are useful in 
understanding the radiative properties of participating species and the nature of nongray 
radiation. The validity of radiative transfer analyses depends upon the accuracy with 
which absorption-emission and scattering characteristics of molecular species are modeled 
and these are reviewed in [12,13]. 

For gray participating mediums between parallel plates and in a cylinder, the exact 
formulations are available in the literature for different physical conditions. Ussikin and 
sparrow [25] studied thermal radiation between parallel plates separated by an absorbing- 
emitting, nonisothermal gray gas. Sparrow, Ussikin, and Hubbard [26] investigated 
radiative transfer in a nonisothermal gray spherical medium. Viskanta [27] studied 
the interaction of conduction, laminar convection, and radiation in a plane layer of a 
radiating fluid. Cess and Tiwari [28] investigated heat transfer to laminar flow of an 
absorbing-emitting gas between parallel plates. Tiwari [29] studied radiative interaction 
in transient energy transfer in gaseous systems. Tiwari and Singh [30] extended Tiwari’s 
work [29] for fully developed laminar flows. The studies presented in [27-29] have 
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reviewed other available literature on gray as well as nongray radiative transfer between 
planar geometries. 

Einstein [31] considered radiant energy transfer in an absorbing-emitting gray gas 
flowing within a black walled cylindrical pipe. Nichols [32] studied the influence of 
the absorption of radiation on the temperature profile and heat transfer to an absorbing 
medium flowing turbulently in an annulus. deSoto and Edwards [33] predicted the 
radiative interchange between a black tube and a nonisothermal nongray gas within the 
tube. deSoto [34] investigated the coupling of radiation with convection and conduction 
in a nonisothermal nongray gas flowing in the entrance region of a black walled tube. 
Pearce and Emery [35] treated the thermal entry region for laminar flow of a gray gas or a 
gas with gray bands. Kesten [36] presented the equation for the spectral radiant heat flux 
distribution in an absorbing-emitting gas contained in a long cylinder. Landram, Greif, 
and Habib [37] studied heat transfer in turbulent pipe flows with optically thin radiation. 
Habib and Greif [38] investigated nongray radiative transport in a cylindrical medium. 
Tiwari and Cess [39] studied heat transfer to laminar flow of nongray gases through a 
circular tube. Wassel and Edwards [40] investigated molecular gas band radiation in 
cylinders for the axial band absorptance, mean beam length, absorptivity and emissivity 
for a nongray gas. 

As mentioned earlier, in case of combustion chambers and industrial furnaces, heat 
transfer takes place by simultaneous convection and radiation. Also, in some postulated 
nuclear reactor loss-of-coolant accident scenarios, it is important to predict heat transfer 
from heated walls to flowing water vapor, which is capable of absorbing and emitting 
thermal radiation. The main emphasis of this study is to present basic formulations and 
solution procedure, to calculate the radiative interaction in an absorbing-emitting medium 
in a laminar flow through parallel plates and a circular duct while other modes of heat 
transfer occur simultaneously. To accomplish this it is essential to review different band 
absorption models available in the literature and incorporate the most suitable model 
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in the governing equations. For certain applications, it is desirable to provide limiting 
forms of the basic governing equations and, if possible, obtain closed form solutions. 
These solutions are very useful in estimating the extent of radiative interaction for a 
given physical problem. 

A review of different absorption models is presented in Chap. 2. The general problem 
of radiative interaction is formulated for flow between black parallel plates in Chap. 3 for 
both gray and nongray mediums. A similar formulation is presented for the flow through 
a black circular duct in Chap. 4. The numerical procedure for both cases is presented 
in Chap. 5. Typical results for both geometries are presented and discussed in Chap. 6. 
Finally, the results and conclusions reached from this study are summarized in Chap. 7. 


Chapter 2 

BAND ABSORPTION AND CORRELATIONS 


The study of radiative transmission in nonhomogeneous gaseous systems requires a 
detailed knowledge of the absorption, emission, and scattering characteristics of the spe- 
cific species under investigation. In absorbing and emitting mediums, an accurate model 
for the spectral absorption coefficient is of vital importance in the correct formulation of 
the radiative flux equations. A systematic representation of the absorption by a gas, in 
the infrared part of spectrum, requires the identification of the major infrared bands and 
the evaluation of the line parameters (line intensity, line half-width, and spacing between 
the lines) of these bands. The line parameters depend upon the temperature, pressure 
and concentration of the absorbing molecules and, in general, these quantities vary con- 
tinuously along a nonhomogeneous path in the medium. In recent years, considerable 
efforts have been expended in obtaining the line parameters and absorption coefficients 
of important atomic and molecular species. 

For an accurate evaluation of the transmittance (or absorptance) of a molecular band, 
a convenient line model is used to represent the variation of the spectral absorption 
coefficient. The line models usually employed are Lorentz, Doppler, and Voight line 
profiles. A complete formulation (and comparison) of the transmittance and absorption 
by these line profiles is given in [9-13]. In a particular band consisting of many lines, 
the absorption coefficient varies very rapidly with frequency. Thus, it becomes very 
difficult and time-consuming task to evaluate the total band absorption over the actual 
band contour by employing an appropriate line profile model. Consequently, several 
approximate band models (narrow as well as wide) have been proposed which represent 
absorption from an actual band with reasonable accuracy [6-13] and [14-24]. Several 


4 



5 


continuous correlations for the total band absorption are available in literature [9-13] and 
[19-23]. These have been employed in many nongray radiative transfer analyses with 
varying degree of success [9-13] and [14-23]. A brief discussion is presented here on 
the total band absorption, band models, and band absorptance correlations. 

The absorption within a narrow spectral interval of a vibration rotation band can quite 
accurately be represented by the so-called “narrow band models.” For a homogeneous 
path, the total absorptance of a narrow band is given by 




J [1 — exp(— fc^X)] du> 


( 2 . 1 ) 


where ku, is the volumetric absorption coefficient, u is the wave number, and X = py 
is the pressure path length. The limits of integration in Eq.(2.1.1) are over the narrow 
band pass considered. The total band absorption of the so-called “wide band models” 
is given by 


OO 

A = J [1 — exp(—k w X)]d(u> — u) 0 ) (2.2) 

— OO 

where the limits of integration are over the entire band pass and u> 0 is the wave number at 
the center of the wide band. In actual radiative transfer analyses, the quantity of frequent 
interest is the derivative of Eqs.(2.1.1) and (2.1.2). 

Four commonly used narrow band models are Elsasser, Statistical, Random Elsasser, 
and Quasi-Random. The application of a model to a particular case depends upon the 
nature of absorbing-emitting molecule. Complete discussions on narrow band models, 
and expressions for transmittance and integrated absorptance are available in the literature 
[9-13] and [14-16]. Detailed discussions on the wide band models are given in [9-13] 
and [17-23]. The relations for total absorptance of a wide band are obtained from 
the absorptance formulations of narrow band models by employing the relations for the 
variation of line intensity as [9-13] and [20-23] 


d 



bp \lo - ui 0 \ \ 

A 0 J 


(2.3) 
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where Sj is the intensity of the jth spectral line, d is the line spacing, S is the integrated 
intensity of a wide band, A 0 is the band width parameter, and b 0 = 2 for a symmetrical 
band and b 0 = 1 for bands with upper and lower wave number heads at u> 0 . The total 
absorptance of an exponential wide band, in turn, may be expressed by 


A(u,0) = _ J_ j [Ajv(u,/3)j d(u> -w 0 ) (2.4) 

wide band 

where u = is the nondimensional path length, /? = is the line structure param- 
eter, 7 l is the Lorentz line half-width, and Apt(u,j3) represents the mean absorptance 
of a narrow band. 


By employing the Elsasser narrow band absorptance relation and Eq.(2.1.3) the 
expression for the exponential wide band absorptance is obtained as [ 12 ] 

T 

A(u,0) = 7 + J [In ip + E x (ip)] dz ( 2 . 5 ) 

o 

where ip = I ) 1 T = 0.5772156 is the Euler’s constant, and Ei(ip) is the 

exponential integral of the first order. Analytic solution of Eq.(2.1.5) can be obtained 
in a series form as [12,13] 


where 


A(uJ) = 

71 = 1 


A n SUM( mn) 
n(B + 1 )" n! (n — 1)1 


( 2 . 6 ) 


SUM(mn) = 

m=0 

A = —u tanh p , B — 


(n + m — 1 )! ( 2 m — \)\C T ‘ 


2 m (m\y 


C = 


cosh / 3 
2 2 B 


(2.7) 


(1 + cosh/?) ( B + 1 ) 

The series in Eq.(2.1.6) converges rapidly. When the weak line approximation for the 
Elsasser model is valid (i.e. j3 is large), then Eq.(2.1.5) reduces to 


A (u) = 7 + In (u) -I- E\ ( u ) 


( 2 . 8 ) 
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In the linear limit, Eqs.(2.1.5) and (2.1.6) reduce to A = u, and in the logarithmic limit 
they reduce to A = 7 + ln(u). It can be shown that Eq.(2.1.5) reduces to the correct 
limiting form in square-root limit. Results of Eqs.(2.1.5) and (2.1.6) are found to be 
identical for all pressures and pathlengths. For p > 1 atm., results of Eqs.(2.1.5)-(2.1.7) 
are in good agreement for all pathlengths. 

By employing the uniform statistical, general statistical, and random Elsasser narrow 
band models absorptance relations and Eq.(2.1.3), three additional expressions for the 
exponential wide band absorptance were obtained in [12,13]. The absorptance results of 
the four wide band models are discussed in detail in [13], The expression obtained by 
employing the uniform statistical model also reduces to the relation given by Eq.(2.1.7) 
for large (3. 

Several continuous correlations for the total absorptance of a wide band, which are 
valid over different values of path length and line structure parameter, are available in 
literature. These are discussed, in detail, in [9-13] and [20-23] and are presented here in 
the sequence that they became available in the literature. Most of these correlations are 
developed to satisfy at least some of the limiting conditions (nonoverlapping line, linear, 
weak line, and strong line approximations and square-root, large pressure, and large path 
length limits) for the total band absorptance [10-13]. Some of the correlations even have 
experimental justification [9-18]. 

The first correlation for the exponential wide band absorptance (a three piece cor- 
relation) was proposed by Edwards et al. [17,18]. The first continuous correlation was 
proposed by Tien and Lowder [9], and this is of the form 


A (u, (3) = In 


uf(i) 


(u 2) 

(« + 2 /( 0 ) 


+ 1 


(2.9) 


where 


/ ( t ) = 2.94 [1 — exp (— 2.60f)] , t = ^ 
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This correlation does not reduce to the correct limiting form in the square-root limit 
[10,13], and its use should be made for /? > 0.1. Another continuous correlation was 
proposed by Goody and Belton [22], and in terms of the present nomenclature, this is 
given by 


A(u,0) = 2 In 



0 = 2 1 


( 2 . 11 ) 


Use of this correlation is restricted to relatively small 0 values [10-13]. Tien and Ling 
[23] have proposed a simple two parameter correlation for A(u,0 ) as 


.4(u) = sinh 1 (u) 


( 2 . 12 ) 


which is valid only for the limit of large 0. A relatively simple continuous correlation 
was introduced by Cess and Tiwari [12,13], and this is of the form 


A(u,0 ) = 2 In 



(2.13) 


where 0 = 41 = By slightly modifying Eq.(2.1.13), another form of the wide band 
absorptance is obtained as [29] 


where 


A (it, 0) = 2 In 



0 . 1 , 0 < 1 and all u values 
c = 0 . 1 , 0 > landu < 1 


(2.14) 


0.25, >3 > landu > 1 . 

Equations. (2. 1.1 1) and (2.1.14) reduce to all the limiting forms [10]. Based on the 
formulations of slab band absorptance, Edwards and Balakrishnan [20] have proposed 
the correlation 


A(u) = In (u) + Ei (u) + 7 + - - E 3 (u) 


(2.16) 
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which is valid for large j3. For present application, this correlation should be modified 
by using the technique discussed in [12,13]. Based upon the formulation of the total 
band absorptance from the general statistical model, Felske and Tien [21] have proposed 
a continuous correlation for A(u,j3 ) as 


A (u, /?) = 2E 1 ( t Pu ) + E X (^) - Ex [(^) (1 + 2 1 


(1 + 2 1 ) 


where 


!+ - 

\u 


(2.17) 


(2.18) 


The absorptance relation given by Eq.(2.1.17) is another simple correlation which is valid 
for all path lengths and for t = (/?/ 2) > 1. The relation of Eq.(2.1.16) can be treated as 
another correlation applicable to gases whose spectral behavior can be described by the 


Elsasser model. In ref. 13, Tiwari has shown that Elsasser as well as random model 


formulations for the total band absorptance reduce to Eq.(2.1.17) for t > 1. 


Band absorptance results of various correlations are compared and discussed in some 
detail in [12,13] and [24], It was found that results of these correlations could be in error 
by as much as 40 % when compared with the exact solutions based on different models. 
Felske and Tien’s correlation was found to give the least error when compared with the 
exact solution based on Elsasser model. The results of Cess and Tiwari’s correlations 


followed the trend of general statistical model. Tiwari and Batki’s correlation [Eq.(21.16) 
or (2.1.17)] was found to provide a uniformly better approximation for the total band 
absorptance at relatively high pressures. The sole motivation in presenting the various 
correlations here is to see if their use in actual radiative processes made any significant 
difference in the final results. 


In Ref. 24, use of several continuous correlations for total band absorptance was 
made to two problems to investigate their influence on the final results of actual radiative 
processes. For the case of radiative transfer in a gas with internal heat source, it 
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was found that actual center-line temperature results obtained by using the different 
correlations follow the same general trend as the results of total band absorptance by 
these correlations. From these results, it may be concluded that use of the Tien and 
Lowder’s correlation should be avoided at lower pressures, but its use is justified (at 
moderate and high pressures) to gases whose spectral behavior can be described by the 
regular Elsasser band model. For all pressures and path length conditions, use of the 
Cess and Tiwari’s correlations could be made to gases with bands of highly overlapping 
lines. In a more realistic problem involving flow of absorbing emitting gas, results of 
different correlations (except the Tien and Lowder’s correlation) differ from each other 
by less than 6 % for all pressures and path lengths. Use of Tien and Lowder’s correlations 
is justified for gases like CO at moderate and high pressures. For gases like CO 2 , use of 
any other correlation is recommended. While Felske and Tien’s correlation is useful for 
all pressures and path lengths to gases having random band structure, Tiwari and Batki’s 
simple correlation could be employed to gases with regular or random band structure 
but for p > 1.0 atm. 



Chapter 3 

LAMINAR FLOW BETWEEN PARALLEL PLATES 

This chapter covers the development of different types of formulation for steady 
laminar flow of absorbing-emitting, constant property, incompressible gas through two 
parallel black plates. Basic formuladon for such a flow is given in Sec. 3.1. Gray gas 
formulation is developed with limiting cases in Sec. 3.2 and Nongray formulation is 

given in Sec. 3.3. 


3.1 Basic Formulation 

The physical model for the parallel plate geometry is shown m (Fig. 3.1.1). 
For radiation participating medium, the equations expressing conservation of mass and 
momentum remain unaltered, while the conservation of energy for such a flow is given 

asfll 

pC p j£- = div (fcVT) + I3T^ + P* - div< l R (3-1) 

where p is the dynamic viscosity, fi is the coefficient of thermal expansion of the fluid 
and <f> is the Rayleigh dissipation function. 

The Condition of uniform surface heat flux for each plate is assumed such that the 
temperature of the plates varies in the axial direction. Fully developed heat transfer is 
considered and axial conduction and radiation is assumed to be negligible as compared 
with the normal components. In other words, this represents physical conditions of a 
large value of Peclet number[l]. Consequently Eq.(3.1.1) can be expressed as 



' dT dT' 

, a 2 T , dp 

du 

pC p 

[ u !h +v -fy\ 

= + + " 

dy 
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If in addition, it is assumed that the Eckert number of the flow is small, then Eq.(3.1.2) 
simplifies to 


dT dT d 2 T 1 dq R 

U dx + V 8y ~ ° dy 2 ~ pC p dy 

where a = (k / pC p) represents the thermal diffusivity of the fluid. 


(3.3) 


For a steady fully developed flow, v = 0 and u is given by the parabolic velocity 
profile as 


u = 6u m (£ - e) ; £ = J 


(3.4) 


where u m represents the mean velocity. Also for the flow of a perfect gas with uniform 
wall heat flux, j s constant and is given by 


dT _ 2aq w 
dx u m LJk 

Now, by combining Eqs. (3. 1.3-3. 1.5), the energy equation is expressed as 

d 2 T _ \2q w , _ 2 \ 1 dqR 

dy 2 L/k ^ ^ > k dy 

Upon defining nondimensional temperature as 


(3.5) 


(3.6) 


e = (T - TQ 

(q w L/k) 


the energy equation, Eq.(3.1.6), is expressed as 


d 2 e 

d( 2 


12 (a -a 2 ) 


1 dqR 
qw d( 


(3.7) 


(3.8) 


The boundary conditions for this problem can be expressed as 

At t he surface ; 0 (0) =#(1) = 0 
, 1 dO n 

AH = 2 ; — = ^ = 0 and q R = 0 

/d0\ _ _ /dd\ 


(3.9a) 

(3.9b) 

(3.9c) 
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It should be noted that all the boundary conditions given in Eqs. (3.1.9) are not 
independent and any two conditions can be used to obtain solutions. 

Equation. (3.1.8) can be integrated once and using the boundary condition, one can 
obtain another form of the energy equation as 

9 ( - 2 (3£ 2 - 2£ 3 ) + 1 = (3.10) 

Qw 


For flow problems, the quantity of primary interest is the bulk temperature of the 
gas, which may be expressed as [11] 

l 

e h = (r * T - P = 6 f 0(0 (£ - £ 2 ) d( (3.11) 

qwL/k J 
o 


The heat transfer q w is given by the expression, q w = h c ( T\ — T b ), where h c is 
the effective heat transfer coefficient (W/cm 2 — K). In general, the heat transfer results 
are expressed in terms of the Nusselt number, Nu = h c D b /k. Here, Dh represents the 
hydraulic diameter, and for the parallel plate geometry it equals twice the plate separation, 
i.e., Dh = 2L. Upon eliminating the heat transfer coefficient h c from expression for q w 
and Nu, a relation between the Nusselt number and the bulk temperature is obtained as 


Nu = 


2Lq h 


k(T x -T b ) 0 b 


(3.12) 


The heat transfer results, therefore, can be expressed either in terms of Nu or #5. In 
order to accomplish this, a proper formulation for the radiative flux appearing in Eqs. 
(3.1.8) and (3.1.10) must be provided. 

As discussed in Chap. 2, the expression for the radiative flux is formulated in 
terms of the absorption coefficient, which in turn is a strong function of frequency. 
Probably the greatest simplification in formulating the radiative flux model is the gray 
medium approximation. This is discussed in the following section. The gray medium 
approximation also serves as an initial step toward nongray analyses. 
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3.2 Gray Gas Formulation 

In gray gas formulation (gray-medium approximation), the absorption and the scat- 
tering coefficients are assumed to be independent of the wavelength. For a gray non- 
scattering medium with black bounding surfaces the equation for radiation flux can be 
given as [1] 


qR = 2aT 4 E 3 (r) - 2aT 2 E 3 ( t - r) 


where 


+ 2a J T 4 (t) E 2 (r -t)dt -2a J T 4 ( t ) E 2 (f - r) dt (3.13) 

0 r 

y L 1 

r = J K(y) dy , r 0 = J n{y) dy , (<) = j y, n ~ 2 e~ » dp 

0 0 0 
In the preceding equations, E n (t) are the exponential integral functions, and r and r G 

represent the optical coordinate and optical path, respectively. For exponential kemal 

approximation, 

E 2 (t) ~ |e-f* ; E 3 (t) ~ (3.14) 

Consequently, Eq. (3.2.13) is written as 


QR 


= aT 4 e~$ r - <7T 2 4 e“5( r °- r ) 


3 _3 

+ 2 ae 2 


• >o 

T J T 4 (t)e^dt - ^ae* r jT 4 (t) 


e * l dt 


(3.15) 


Upon differentiating Eq. (3.15) twice, the integrals can be eliminated and the 
following expression is obtained 


d 2 q R 9 0 dT 4 

= 3,7 — 

For linerized radiation Eq.(3.16) becomes 


jqR 

dr 2 


9 o 

-qR = YlaTl 


dT 

dr 


(3.16) 


(3.17) 
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where k p is the Planck mean absorption coefficient, Eq. (3.17) can now be written as 

_ £ t 2 ? = (3.19) 

dt? 4 oQR [ kK p J q d( 

In obtaining Eq.(3.19), it has been assumed that the absorption coefficient (not the Planck 
mean absorption coefficient) is not a function of temperature. This is consistent with the 
assumption of constant fluid properties in the medium. The boundary conditions for 
Eq.(3.19) are found to be 



(dqR) 

V d( )t =0 


(3.20) 


Using the dimensionless quantities defined by Eq. (3.18) the energy equation given 
by Eq. (3.8) is expressed as 

(3.21) 

? ' di 2 q di 

For this case, the boundary conditions become 


0(0) = 0 




df, \2 


(3.22) 


Equations (3.19) and (3.21) can now be solved simultaneously and the following 
result is obtained for the bulk temperature 

0 b = Ci [24 - 12 M\ + Ml + [Ml - 12 M x - 24) e _Ml l 
__ 12ji_ 1771 _ 17 

5A/j 4 + 70 A/ 2 70 


(3.23) 
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where 


Ml = 

4 » 

+ 71 


12<rT 3 r 2 

W O 

71 = 

k « 

P 

n, — 

71 


L-l — 

Mf 

>o(l 


48 — 3r 0 A/j 4- 36 r 0 
- e~ Ml ) + 2Mi (1 + e~ Ml ) 


(3.24) 


The bulk temperature can be easily calculated from Eq.(3.23) for any gas once the 
Planck’s mean absorption coefficient k p is known. 

Two limiting solutions can also be found for large and small optical thicknesses. 
These are called optically thin limit and optically thick limit solutions. 


3.2.1 Optically Thin Limit 


The solution in the optically thin limit is obtained by making use of the parameter 
Mi defined in Eq. (3.24) as 

Ml = - a t 2 0 + 71 (3.25) 

In the optically thin limit r 0 — > 0 and, therefore. 


lim (Ml) - > 7i , lim (Mi) = 

> 0 r G — >0 


(3.26) 


Consequently, 


24 

lim (Ci) = — 

r ° 7i (1 + e~^) 


(3.27) 


The expression for the bulk temperature in thin limit, therefore, is obtained as. 


i^b)thm l 


24 


24 


(f) , 1 

( 7? I s/li 


7i (l + e 

1 _ e -v^ 


- 12 VtT+ (7? - 12^-24) e-v^ 
! + e -V7TJ - 288 + 24 ^ - 2 ^} 


12 

571 


(3.28) 
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3.2.2 Optically Thick Limit 


In the optically thick limit r 0 — > oo. From the definitions of 7 1 and Mi given 
in Eq. (3.24), it is obvious therefore t 0 >> 1, M\ >> 1. Equation (3.23), therefore, 
reduces to 


__H+I I 

( b)thxck 7Q + 7Q 


iY(3 + » 

and rearranging the terms, a simpler form is obtained as 

1 


thick 7 q 


1 + (3^) 


where, 




4<J T* 


(3.29) 


(3.30) 


3.3 Nongray Formulation 

In this section, formulation of the problem under realistic conditions is tried. Under 
realistic conditions the gas is not gray, that is the band parameters are functions of 
wavelength. Now, for solving Eq.(3.1.10), first we need to get equation for radiative 
flux. The radiative transfer equations are formulated for one-dimensional planar systems. 
For diffuse boundaries and in the absence of scattering, the expression for the radiative 
flux is given as [1] 


Qrx (t x ) = 2B 1X E 3 ( t a ) - 2 B 2 \E $ (t oX - t x ) 


+ 2 


where 


J e 6 A (0 E 2 (t\ — t)dt — J e bX ( t ) E 2 ( t-r \ ) 

-0 Tx 

y L 

T \ = J k x dy ■ t oX — J k x dy 

0 0 

1 

E„ (t) = J fi n ~ 2 e~idfi 


dt 


(3.31) 


(3.32a) 


(3.32b) 
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In the preceding equations, E n (t) are the exponential functions, and r A and r oA 
represent the optical coordinate and optical path, respectively. The quantities B 1A and 
B 2A represent the spectral surface radiosities, and for nonreflecting boundaries, Eq.(3.31) 
is expressed in terms of the wave number as [35] 




+ 2 


• (Jj 1 OUJ 

J Flw (*) E 2 (t u — t)dt — j F 2w ( t ) E 2 ( t - T w ) dt 


(3.33) 


where 


Flw (0 (0 j F (t) (0 ^2w 

Equation (3.33) is the general equation for one-dimensional absorbing-emitting 
medium with diffuse non-reflecting boundaries. For nongray analyses, it is often con- 
venient to replace the exponential integrals by appropriate exponential functions. Upon 
employing the exponential kemal approximation 


E 2 (t) - fexp(-|i) , E x (t) = §exp(-§<) 
Equation (3.33) is expressed in physical coordinates as 


?Pw(y) = - e 2ui 


+ 


if 

^ J F Xu ,(z)k u exp 
0 

L 

^ J F 2w ( 2 ) k u exp 


- z) 


~2 kuJ ( Z ~ y ) 


dz 


dz 


(3.34) 


y 

where z is a dummy variable for y. In obtaining Eq. (3.34) it has been assumed that 
temperature differences in the medium are small and k u is independent of temperature. 
This is consistent with the assumption of constant properties within the medium. 
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As discussed in Chap. 2, the total band absorptance can be expressed as 

OO 

A(y ) = J [1 — exp(—k u y)]diL> ~ cm -1 (3.35) 

o 

where both k w and u> have units of cm -1 . Differentiation of Eq.(3.35) gives 

OO 

A 1 (y) = J k^exp^—k^y) du> ~ cm~ 2 (3.36) 


and 


CC 

A" (y) = J -klexp(-k u y)du 


~ cm 


-3 


(3.37) 


Equations (3.35)— (3.37) are employed to express Eq.(3.34) in terms of the band absorp- 
tance. 

The total radiative flux is given by 


such that 


OO 

?*(y) - J qRu,(y)du 


(3.38) 


dqR(y) f dq Rt 


dy 


f dqRw d f 

“ J - Ty J 


(3.39) 


o o 

Upon substituting Eq.(3.34) in Eq.(3.38) there is obtained for a multiband gaseous system 


QR{y) =e l-e 2 + ^]T f 

| J F iu , (z) k u ,exp ( y - z) 

^ 0 

L 

- J F 2 u>, (2) k Ut exp 


(3.40) 


~ 2 kw ' ( z ~ y) 


dz 


dz \duj % 


It should be noted that the following relations have been used in obtaining Eq.( 3 . 40 ) 


OO OO 

fe lui dcj = ej ; f e 2u> du> = e 2 
0 0 
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OO 


I l f F lu (z)k u exp (y - z) 

0 l 0 

E/ 1 / ( z ) k *. ex P 


dz 


} duj = 


~2 k “' ( y ~ z ^ dz 


> dui 


‘Au, VO 

where n represents the number of bands in a multiband system. 

By utilizing the definition of the band absorptance and its derivatives as given by 
Eqs.(3.35)-(3.37) and evaluating the value of the Planck function at the center of each 
band, Eq.(3.40) is expressed as 

y 


qR(y) = ei - e 2 + 2 E i /*u.(*)4 


1=1 


!<»-*> 


L 

~ J Fv Wot (z) A, 

y 

where u> 0 i represents the center of the ith band. 

Equation (3.41) is in proper form for obtaining the nongray solutions of molecular 
species. However, in order to be able to use the band model correlations, this equation 
must be transformed in terms of the correlation quantities defined in Chap. 2. The 
following quantities, therefore, are needed for the transformation 




dz 


dz 


(3.41) 


u — 


(l) ra; ""=(l) p£;p5= /^ 


(3.42) 


Now, by using the definition A = - j -, Eq.(3.36)-(3.37) are written as, 
ji . . d rj , 1 dA(y) A' (y) , 


A'(y) = A^=A, f^l^ . 

dy [ du dy ] 

A"(y) = [PS(T)] 2 (F)j' 

\ ** O / 


= PS(T)A'(u ) 


(3.43a) 

(3.43b) 

(3.43c) 


The dimensions of both sides in Eq.(3.43b) and ( 3.43c) agree with the dimensions 
given in Eqs.(3.36) and (3.37). By employing the definitions of Eqs.(3.42) and (3.43), 
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Eq.(3.41) is expressed as 


n ( u> r 

(«) = e l _ e 2 + | YL A °*\f ( U «) A'i | («. - «<) 

1=1 ^ o 

«ot \ 

- J F 2Ui (u[) A[ | (u' - tii) du\ | (3.44) 


where u' is the dummy variable for u and A' (it) = j£. It should be noted that F lwi and 
F^j in Eq.(3.44) represent the values of F lu , and F 2w at the center of the ith band, and 

3? = (d?) (ft) = [tP] {%)■ 


By defining the new independent variables as 

, - — - L . A - — - 1 

S r 5 S T 

u 0 L u 0 L 

Eq.(3.44) can be expressed as 


(3.45) 


QR (0 = e i - e 2 + | ^2 ^o.Uoij J F lui (£') A\ |t i 0l (Z - Z') dZ' 

»=i l o 

- J F 2u , (?) % d? | (3.46) 


where again A' (u) denotes the derivative of A(u) with respect to u, and = 

(t) (&) = (i) (t)- 


Equations (3.44) and (3.46) allow us to make use of the band model correlations for 

the wide-band absorptance because these correlations are expressed in terms of u and j3. 

However, it is often desirable and convenient to express the relations for qR and div qR 

which only involve A (u) and A' ( u ) but not A" ( u ). This is accomplished by integrating 

the integrals in the above equation by parts. This results in simpler integrals. Upon 

performing the integration, Eq. (3.46) is expressed as 

, i . 


9R( 0 = e l 


- e2+ ip m {j[±^P]A.[^ i{( -?)\,? 

- / Ai [I" 1 ” “ () ] <3 ' 47) 
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By differentiating Eq. (3.47) using the Leibinitz formula, there is obtained 

L 1= 

dt 


1 = 1 


de u , (?) 


i 

/[ 


d? 


de u , (£') 


A' 


2 u <» (f ~ O 


d(' 


A'. 


^ ot (C - o 


d? 

d? 


(3.48) 


Equations (3.47) and (3.48) are the most convenient equations to use when employing 
the band-model correlations in radiative transfer analyses. 

For the present physical problem, ei = e 2 and F lu;i = F 2wi . Thus, for the case of 
linearized radiation, a combination of Eqs.(3.10) and (3.46) results in [29] 


O' - 2 (3£ 2 - 2£ 3 ) + 1 » | (j) 

\ u oi {( - O 


1 

- f no % 

£ 








(3.49) 


where H u = A oi (T) ; i?i = £ H u 

«=i 

A combination of Eqs.(3.10) and (3.48) gives an alternate form of the energy equation 
for the steady case as 


O' - 2(3£ 2 -2£ 3 ) + 1= 



(3.50) 
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Note that this equation can also be obtained directly by integrating the left hand side 
of Eq. (3.49) by parts. Equations (3.49) and (3.50) provide two forms of the energy 
equation for the steady-state conditions. 

For the case of negligible radiation, Eqs. (3.49) and (3.50) reduce to 

O' = 2 (3£ 2 - 2£ 3 ) - 1 (3.51) 

The solution of Eq.(3.51) is found to be 

id) = {(2« 2 -{ 3 -l) (3-52) 

Thus, a combination of Eqs. (3.11) and (3.52) gives the result for the bulk temperature 
for the steady case with no radiation as 

9 b = ~ = -0.24286 (3.53) 

This result is useful in determining the extent of radiative contributions. 


3.3.1 Optically Thin Limit 


In the optically thin limit [1], A(u ) = u, -A'(u) = 1, and A" ( u ) = 0. In this limit, 
therefore, Eqs. (3.49) and (3.50), reduce to 

* i 


O' - 2 (3£ 2 — 2£ 3 ) + 1 = j-N 


J e («') d('-Je (O JC 


L0 


(3-54) 


O' - 2 (3£ 2 — 2£ 3 ) + 1 = -N 




(3. 


55) 


where N = (nr) A'i - (nr) E s i ( T ) (%^) 

Hu = ^«(T)(>) ri ; H i = ±H U 


T, 
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The differentiation of Eqs.(3.54) and (3.55) yields the same energy equation for the 
optically thin limit as 


0" - 3 N6 = 


(3.56) 


The solution of Eq.(3.56) satisfying the boundary conditions 0(0) = 0 and 0(1) = 0 is 
found to be 


*«) 


= (- 


16 

smh ( 

3 iV 2 ) 

sinh 

(*)< 



cosh 




(3.57) 


Alternatively, the solution of Eq.(3.57) is written as 


0(0 = C\ exp + C2exp (— \/rn() 

+ (24 — 12 m£ + 12 m£ 2 ) ; m = 3 N 


(3.58) 


The constants Ci and Cj arc obtained by using the boundary conditions 0(0) = 0 and 
6' (1/2) = 0, and the solution for 0(£) is found to be 


*(0 = 


/ 1 \ f ( 22-=^ f e-v^ e ^ + e -S%) 

H V 1 + e~^J \ t 

\ m J o 

' ' I i a -i d _ c , in . /-z 


(3.59) 


4 " 24 — 12 m£ -f" 12 m £ 2 
Equations (3.59) and (3.3.60) should produce identical results. The expression for the 
bulk temperature, in this case, is obtained by combining Eqs.(3.11) and (3.59) as 


576 

z 

ms 


(3 N) 


1 - e~v^" 

1 + e~v^ 

f 576 

\Vsn 


288 24 12 

m 3 m 2 5 m 


[ 4 =NEXP - 21 . 6 iV 2 + 72 N - 288 1 


(3.60a) 

(3.60b) 


where 


N EX P = 
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3.3.2 Large Path Length Limit 

In the large path length limit, u 0 j » 1 for each band of interest and A(u) = 
In (w) , A' (u) = fa and A" (u) = —fa [29]. Thus, in this limit Eqs. (3.49) and 
(3.50) reduce to 

i 

&' -2 (3£ 2 - 2£ 3 ) + 1 = M J (3 ' 61) 

0 


O' - 2(3£ 2 -2£ 3 ) + 1 = (f)E Fl * 



£ 


where M = Zjfa = (f ) £ A 0 , (fy) Ti 

1=1 

Through integration by parts, it can be shown that Eq.(3.62) reduces to Eq.(3.61). The 
non-dimensional parameter M constitutes the radiation-conduction interaction parameter 
for the large path length limit. Equation (3.61) does not appear to possess a closed form 
solution; a numerical solution, however, can be obtained easily. 



Chapter 4 

LAMINAR FLOW THROUGH CIRCULAR DUCT 

This chapter covers the development of different types of formulation for steady 
laminar fully developed flow of absorbing-emitting, constant property, incompressible 
gas in a circular duct. Basic formulation for such a flow is given in Sec. 4.1. Gray gas 
formulation is developed in Sec. 4.2 and nongray formulation is given in Sec. 4.3. 


4.1 Basic Formulation 


The physical model for the circular duct geometry is shown in (Fig.4. 1). For radiation 
participating medium, the equations expressing conservation of mass and momentum 
remain unaltered, while conservation of energy for such a flow is given by the same 
expression as Eq.(3.1), i.e., 

P C v~^ = div(kVT) + ~ divqn (4.1) 

Under similar conditions as mentioned in Sec. 3.1, Eq.(4.1) can be written in circular 
coordinates as 


dT a d f dT\ 115. , 

U dx ~ r dr V dr ) ~ pC p r dr 

By noting that for a uniform wall heat flux and fully developed heat transfer 

dT _ 2a gu, 
dx u m r 0 k 

and by employing the parabolic laminar velocity profile for u given by 



(4.2) 


(4.3) 


(4.4) 
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Eq. (4.2) can be written as 


r\ k d ( dT\ Id . , iA - 

— = — 7 T r— — (rq R ) (4.5) 

r 0 ) q w dr \ dr J q w dr 

By noticing that for the case of a circular duct u = C„Pr and u 0 = CgPr 0 , Eq. (4.5) 


1 - 1 - 

r„ 


is expressed as 



Defining non-dimensional quantities 


(4.6) 


r u 

r 0 u 0 


9 = 


T-r w 

q w r 0 /k 


kCjP 

7 


(T - T w ) 


Eq.(4.6) can be rewritten as 

The boundary conditions for this case are given as 


(4.7) 


(4.8) 


6 ’ ( 0 ) = 0 

0(1) =0 (4.9) 


Integrating Eq.(4.8) once and by noting that at £ = 0, (d8/d() = 0 and q R = 0, 
we get 


0t + e-2(; = 


gR (0 

Qw 


For the circular geometry, the bulk temperature is defined as 

2t r 0 


6t ~ u^rrj / /” (r > < ( r>r "* 


0 0 

A combination of Eqs. (4.4) and (4.11a) results in 

i 


*6 = 4 J 


(4.10) 


(4.11) 


(4.12) 


The Nusselt number based on the hydraulic diameter Dj, = 2r 0 is given as 

2r 0 q w 2 
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4.2 Gray Gas Formulation 

For a gray non-scattering medium within a black bounding surface, the equation for 
radiative flux can be written using the exponential approximation as [29] 


d_ 

dr 


r dr 


9 2 „ dT 4 

- ^ p qR = 3(7 K p~^t 


For linearized radiation Eq.(4.13) can be written in dimensionless form as 


d 

da 




9 2 dO 

~jr 0 q R = l2 q w — 


where 


72 = 


1 2 ( 7 KpT^rl _ \2aTirl 




> 'o 


Tfy K pr 0 


Boundary conditions for Eq.(4.14) are 


, , , 3 1 

At £ = 1, -q R = 


2 111 t 0 [idi 

q R { o) = o 


A (f «* ) 


(4.14) 


(4.15) 


(4.16) 


(4.17) 


Equation (4.8) can be integrated once to give 

dO qR(0 




Qw C 


Equation (4.14) and (4.16) are combined to give 

(Mk 2 + 1) « = (2« 1 - « 5 - c,() 

where 

M 2 = \ T ° + 72 

The solution of Eq.(4.18) is found to be 

q R = C 2 Ii (M 2 Q + CzKi (M 2 0 + ^ ( MU 3 + - 2MU) 

For a finite solution as £ — > 0, C 3 must be taken to be zero and 

r2< 


qR = C2I1 (M2O + ( M 2 V + H - ™U) 


(4.18) 


(4.19) 


(4.20) 


( 4 . 21 ) 
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From Eqs.(4.16) and (4.20), there is obtained 


c,h (M 2 ) + (8 - Ml) 




I\ (M 2 ) + (“ii (Af 2 o) _ 1 


+ 


16-/2^ 

m 4 


(4.22) 


or 


3r 0 C 2 /i (M 2 ) + (24r 0 - 3 t 0 M 2 2 ) 


= — 2C 2 


/j(m 2 ) + (m 2 j 0 (m 2 o- 


/iW \ 

* A-i. 


3272 


M 4 


$4.23) 


C 2 [3r 0 /j (Af 2 ) + 2 M 2 I 0 (M 2 )] = [3r 0 M 2 2 - 24r 0 - 32] (4.24) 


/ 72?u>\ 

Vm 2 V 

From Eq.(4.2.17) and (4.2.20) one finds 


C 2 = 


3 r 0 M| - 24 r 0 - 32 
2M 2 / 0 (M 2 ) + 3 T 0 h (M 2 ) 


^ = ——I\ (M 2 () + -2}; (Mfc 3 + 8£- 2 MU) +2 e - e 3 - y- (4.25) 

where from the boundary condition 9' (0) = 0, Ci = 0 
Equation (4.24) is integrated once to give 


6 - 


C 2 

q w M 2 


i 0 (m 2 0 + 


72 




+ 4 e-Mio+e, 


ie) 


- - <?. 


(4.26) 


From the boundary condition 0(1) = 0, one finds 


^ c 2 r . 72 ^16 — 3 M*\ 

Ci = ^ (m)+ wA—r^) 


Consequently, Eq.(4.25) becomes 



q w M 2 

+ (Mi? + 16« 2 - 4 Mi? + 3 Mi - 16) + ? - - \ (4-27) 
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The expression for the bulk temperature can be easily obtained as 

I6J1 (M 2 ) 


0 b = C 


8 - M\ 

M\ 


)u 


+ 


(M 2 ) 

11 


Ml 


11 7 2 8 72 

24 M| 3Mj 24 


(4.28) 


where 


C = 


72 


3t q MI — 24 r 0 — 32 


M| L2M2/,, (M 2 ) + 3r 0 J! (M 2 ) J 


Ml = 4- T2 


72 = 


\2aK p Tlr 2 0 


12<r7*r, 2 

Ar/Cp 


t 0 = «»r, 


p' 0 


4.2. 1 Optically Thin Limit 


Recall the definition of M 2 as 

ry 9 O 

M] = jr 2 + -n 

For optically thin conditions r 0 — > 0 and, therefore. 


limr 0 — > 0 (Ml) = 72 , lim r 0 — > 0(M 2 ) — 7| 


limr 0 — > 0(C) — — 


32 

Zyfnlo {yfn) 


16 

yfnh (\fn) 


The expression for the bulk temperature in the optically thin limit, therefore, becomes 


(&t>\hin ~ 3 

7 2 


256 h (y^) 
\/72 lo (yfii) 


- 128 + I672 - 



(4.29) 


By employing the series expansion of the modified Bessel functions for small values of 
the argument 72, Eq. (4.29) can be written in an alternate form, and by letting 72— > 0 
in that form there is obtained the result in the transparent limit (0j)<Ai'n = — 
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4.2.2 Optically Thick Limit 


The asymptotic series for modified Bessel’s functions for large value of M 2 are 

Mi 


Io (M 2 ) ~ 


h (M 2 ) ~ 


\j2i TM 2 
e M ? 

\/2nM2 


1 9 

1 + rrr + + ••• 


1 - 


8M 2 128A/| 

15 


+ ... 


(4.30a) 

(4.30b) 


8 M 2 128 Ml 

Substituting for Io(M 2 ) and Ii(M 2 ) from Eq. (4.30) into Eq. (4.27) and taking the limit 
as t 0 — > a, we get 


(Qb)th,ck ~ 


11 72 8 72 11 

24 ~ 3MJ - 24 


(4.31) 


Defining 


N = kKp 


4*T3 


72 = ^ and Mi = ? (3 + 1) r» 
Eq. (4.31) is expressed in an alternate form as 


thick 


11 

24 


1 + (w) 


(4.32) 


Equation (4.32) also reduces to the correct form for the case of no radiative interaction. 


4.3 Nongray Formulation 

In this section, formulation of the problem under realistic conditions is presented. 
Under such conditions the gas is nongray, that is the band parameters are functions of 
the wavelength. For solution of Eq.(4.10), an appropriate equation for radiative flux is 
needed. For a circular geometry, the expression for the radiative heat flux in terms of 



34 


non-dimensional path length is given as [39] 

1 L 



cos 7 


( u — u')j 

- J [Buc ( u ') - B wc (T^)] A' { ~^ (u + u' — 2 u sin 7 )^ 

U 

+ A ' 


(4.33) 


where a = 1 and b = 5/4. Equation (4.33) is expressed in an alternate form as 
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In order to combine Eq.(4.10) and (4.34), the Planck’s function is linearized as. 
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Now, a combination of Eqs.(4.10), (4.34) and (4.35) yields 


(4.35) 


d9 

di 



d£ 


£sin 7 


35 


i 

- 1 no & («*-{) 


cos 7 




+ j 0 (£')Ai + sin 7 ) df'J d ~1 ( 4 - 36 ) 


i sin 7 

For nongray gases with n-bands, Eq.(4.36) is written in the following form, 
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where 


H ' = A ot {~ 3 f-) Tw ’ A i = Xl 

Combined solutions of Eqs. (4.37) and (4.11) are obtained by numerical procedures. 


4.3.1 Large Path Length Limit 


As noted earlier, in the large path length limit Uoi » 1 for each band of interest, and 
A* (u,) = 1 / tij. In this limit, Eq. (4.37) reduces to 
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The solution of Eq.(4.38) or (4.39) can be obtained either numerically or in closed form. 


Chapter 5 

METHOD OF SOLUTION 


The solution procedure for both cases, the parallel plate geometry and the circular 
duct, are presented in this chapter. In principle, the same numerical procedure applies to 
both the general and large path length limit cases for both geometries. 


5.1 Parallel Plates 

The general solution of Eq.(3.49) is obtained numerically by employing the method 
of variation of parameters. For this, a polynomial form for 0(Q is assumed in powers 
of £ as 

0(O=it a ™Z m (5-1) 

m= 0 

By considering a five term series solution (a quartic solution in £) and satisfying the 
boundary conditions 8 (0) = O' ( 1/2) = 0 and O' (0) = —O' (1), one obtains 

0 (£) = G1 (£ - 2£ 3 + £ 4 ) + a 2 (£ 2 - 2£ 3 + £ 4 ) (5.2) 


and 


O' (£) = fll (1 - 6£ 2 + 4£ 3 ) + a 2 (2£ - 6£ 2 + 4£ 3 ) 
A substitution of Eq.(5.3) in Eq.(3.49) results in 
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(5.4) 


where expressions for 0(£') are obtained from Eq.(5.2) 

The two unknown constants aj and a 2 in Eq.(5.4) are evaluated by satisfying the 
integral equation at two convenient locations (£ = 0 and £ = 1/4 in the present case). 
The entire procedure for obtaining ai and a 2 is described in [35]. With known values of 
ai and a 2 , Eq.(5.2) provides the general solution for B((). The expression for the bulk 
temperature is obtained by combining Eqs.(3.11) and (5.2) as 


Ob 


\la\ -f- 3o 2 
70 


(5.5) 


The governing equation for the large path length limit is Eq.(3.61). For this equation 
also the solution is given by Eqs.(5.2) and (5.5) but the values of a’s are different in 
this case [35], 


5.2 Circular Duct 

The general solution of Eq.(4.32) is obtained numerically by employing the method 
of variation of parameter similar to that discussed in Sec. 5.1. For this, a polynomial 
form for 0(£) is assumed in powers of £ as given by Eq. (5.1). By considering a 
five term series solution (a quartic solution in £) and satisfying the boundary conditions 
0(1) = 0, 0' (0) = 0 and 0(1) = 0(— 1), one obtains 

0(0 = ai(0-l)+a 2 (0-l) (5.6) 


and 


0'(O = 2 ai £ + 4a 2 £ 3 


(5.7) 
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A substitution of Eqs.(5.6) and (5.7) in Eq.(4.32) results in 
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where expression for 6 (£') are obtained from Eq.(5.7). 

The two unknown constants ai and a 2 in Eq.(5.8) are obtained by satisfying the 
integral equation at any two locations of £. The entire procedure is described in Appendix 
A and relations for the constants are expressed as 


a\ — ((22104) — (ZZ 2 Q 2 )) /DEN (5.9) 


a 2 = {{zz 2 a x )-{zz x a z ))/DEN (5.10) 


DEN = (aic*4 — 0203) 


(5.11) 


The quantity appearing in Eqs.(5.9)-(5.11) are defined in Appendix A. 

Now, with known values of ai and a 2 , Eq.(5.6) provides the general solution for 6(0- 
The expression for the bulk temperature is obtained by combining Eq.(4.1 1) and (5.6) as 



(5.12) 
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The governing equation for the large path length limit is Eq.(4.34). For this equation 
also the solution is given by Eq.(5.6) but the values of ai and a2 are different. In large path 
length limit, the integrals can be evaluated in close forms. Procedure is described briefly 
Equation (4.34) is expressed in an alternate form as 
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Substituting for #(£') in terms of ai and a2 as given by Eq. (5.6), and rearranging 
Eq.(5.14), one obtains 
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By evaluating the inner integral Eq.(5.15) is expressed as 
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Now solving these trigonometric integrals and simplifying, one obtains 
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Equation (5.2.17) can now be written as 
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At £ = 1, a similar solution procedure yields the result 
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Now, with known values of ai and a 2 , Eq. (5.6) provides the general solution for 
0(£) and the bulk temperature is given by Eq. (5.12). 



Chapter 6 

RESULTS AND DISCUSSION 


Extensive results have been obtained for variation of temperature ( 9 ) and bulk 
temperature (#b) at different conditions for all four species in both the cases, and most 
of these are available in Appendix B and C. The computer programs used for numerical 
solution for both geometries are provided in Appendix D and E. Selected results are 
presented here to compare and illustrate the variation of temperature 9 with £ and of 
bulk temperature 9 \> with plate spacing L or duct radius r 0 . Variation of bulk temperature 
with pressure is also presented to illustrate the effect of large pressure path length and 
it’s approach to limit of large Uoi much more clearly. 

Although the results are similar in nature for both the geometries, they are presented 
in seperate sections for clarity. 


6.1 Parallel Plates 

The results for temperature variations across the duct are presented in Figs. 6. 1-6.4. 
Figure 6. 1 shows temperature profile for H 2 O at T w = 500 K, P = 1 atm and L = 10 cm 
with general as well as in the limit of large u 0 j. This demonstrates that the limit of large 
u 0 j overestimates the radiative energy transfer giving a conservative approximation. Thus, 
for practical and realistic problems, it might be easier to find out the extent of radiative 
interaction using the limit of large u Q j formulation and then use the general formulation 
which is computationally complex and expensive. Figure 6.2 demonstrates the effect 
of higher wall temperatures on temperature profiles. As the wall temperature increases 
temperature profile becomes less and less parabolic. In other words, temperature gets 
distributed more and more evenly across the duct. The effect of increasing pressure at any 
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temperature shows a similar trend (Fig. 6.3). At high pressures, the results approach the 
limit of large u<,j. The comparison of temperature profiles for the four species presented 
in Fig. 6.4 demonstrates the relative importance of each species in energy transfer. Water 
vapor (H 2 O) is the most radiation participating species among the four considered here. 
For the physical conditions of Fig. 6.4, OH is the least radiation participating species. 

The bulk temperature results as a function of the distance between the plates are 
presented in Figs. 6.5-6.9. The limiting value of 9 \, = —0.243 corresponds to negligible 
radiation. For a particular plate spacing L, the large path length results represent the 
limiting solutions for high pressures. The results, in general, demonstrate that the effect 
of radiation increases with increasing plate spacing. 

General as well as limit of large u 0 i solutions for the bulk temperature are illustrated 
in Figs. 6.5-6.7 for H 2 O. The results for H 2 O at T w = 1000 K are presented for P = 1, 
2, 5, 10 atm. It is noted that as pressure increases the ability of gas to transfer radiative 
energy increases, approaching the correct limiting solution of large Uoi. The results shown 
in Fig. 6.6 for H 2 O, demonstrate the effect of increasing wall temperature. It is obvious 
that the radiative transfer is more pronounced at the higher wall temperature. The results 
presented in Fig. 6.7 illustrate both effects. 

Comparison of bulk temperature results with plate spacing for various species is 
presented in Figs. 6.8 and 6.9 for a pressure of one atmosphere. The results clearly 
demonstrate the relative ability of the four species for radiative transfer at different path 
lengths. For lower plate spacings and relatively higher temperatures, however, CO 2 
shows a significantly higher ability than other species. This is a typical distinguishing 
feature of the CO 2 under optically thin conditions[29]. 

The trend of the general solution approaching the limit of large u c i (large path length 
limit) is illustrated in Fig. 6.10. A comparison of results for various species is shown 
in Fig. 6.11. 
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Figure 6.5 Variation of bulk temperature with plate spacing for H 2 O; 


T w = 1000 K. 
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Figure 6.7 Variation of bulk temperature with plate spacing for H 2 O. 
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6.2 Circular Duct 

All results for circular duct are presented in Figs. 6.12-6.22. The results are for the 
same species at the same conditions as for the parallel plate geometry. The results, in 
general, exhibit the same trend as presented in Figs. 6.1-6.12 for the parallel plate 
geometry. However, the extent of radiative interactions is entirely different. This 
is because the circular geometry provides additional degrees of freedom for radiative 
interactions[39]. Other basic differences are noted below. 

Results of temperature variations across the duct are plotted from £ = -1.0 to 1.0, 
because of the difference in the location of the coordinate axes as shown in Fig. 4.1. The 
limiting value of #b = -0.4583 corresponds to the case of negligible radiation. Another 
important point to note is that the general solution approaches the limit of large Uoi at 
lower path lengths than for the parallel plate geometry. Extensive results for circular 
duct for all species were obtained. Some of these results are included in Appendix C. 
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Figure 6.14 Comparison of temperature variations across the duct for H 2 O; 


T w = 1000 K and r 0 = 5 cm. 






Figure 6.18 Variation of bulk temperature with radius for H 2 O. 













Chapter 7 

SUMMARY AND CONCLUSIONS 


A brief review is presented on various band models and band model correlations 
that are useful in nongray radiative transfer analyses. Different formulations for one- 
dimensional radiative flux are provided. These are used to develop the basic governing 
equations for energy transfer in gaseous systems. 

Analytical formulations and numerical procedures have been developed to investigate 
the radiative interaction of absorbing-emitting species in laminar fully developed flows 
between parallel plates and through a circular duct. Extensive results have been obtained 
for OH, CO, CO 2 , and H 2 O for different physical conditions. Illustrative results for the 
temperature distribution and bulk temperature arc presented for different pressures and 
wall temperatures. The general nongray results for the circular duct have been obtained 
for the first time. In these results, a lower value of temperature implies a higher ability 
of the gas to transfer radiative energy. 

The gray and nongray formulations for circular duct has been derived in local 
thermodynamic equilibrium. This formulation involves four integrals of which two are 
angular, one is spatial, and one is spectral. One angular integral is taken care of by using 
exponential kemal approximation [27, 32, 37]. Furthermore, the spectral integration for 
nongray formulation is represented by using a total band absorptance. The remaining two 
integrals, one angular and one spatial, are computed numerically by multidimensional 
quadrature method using Gauss-Kronrod rule [41]. The resulting equations are either 
integral or integro-differential and cannot be solved analytically. They are therefore 
solved numerically using the method of variation of parameter [1, 35]. 
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The results, in general, demonstrate that the effect of radiation increases with 
increasing plate spacing/radius, and the radiative transfer is more pronounced at higher 
wall temperature and pressure. One important fact is that, for all species, the general 
solution approaches the correct limiting solution for large u G j. Of all the species, HiO is 
a highly radiation participating species (as compared to C0 2 , CO and OH). The extent 
of radiative interactions is higher in circular duct, because it provides additional degrees 
of freedom. 

The present study provides different kinds of limiting solutions for both the geome- 
tries. These limiting solutions are in closed form and therefore computationally less 
expensive to obtain. In many practical and realistic problems the limiting solutions are 
very useful to demonstrate the importance of radiative interaction. Finally it is important 
to note that any kind of extensive radiation model (line-by-line or narrow band model) 
can be easily used in the nongray formulation provided in this study. 



69 


REFERENCES 

1. Sparrow, E.M. and Cess, R.D., Radiantion Heat Transfer , Brooks/Cole, Belmont, 
Calif., 1966 and 1970. New Augmented Edition, Hemisphere Publishing Corp., 
Washington, D.C., 1978. 

2. Hottel, H.C. and Sarofim, A.F., Radiantive Transfer , McGraw-Hill Book Co., New 
York, 1967. 

3. Siegel, R. and Howell, J.R., Thermal Radiantion Heat Transfer , Mcgraw-Hill Book 
Co., New York, 1971; Second Edition, 1981. 

4. Ozisik, M.N., Radiantive Transfer and Interaction with Conduction and Convection , 
John Wiley & Sons, Inc., New York, 1973. 

5. Edwards, D.K., Radiantion Heat Transfer Notes , Hemisphere Publishing Corporation, 
Washington, D.C., 1981. 

6. Cess, R.D., “The Interaction of Thermal Radiantion with Conduction and Convection 
Heat Transfer”, Advances in Heat Transfer , Vol. 1, Academic Press, New York, 1964. 

7. Sparrow, E.M., “ Radiantion Heat Transfer between Surfaces”, Advances in Heat 
Transfer , Vol. 2, Academic Press, New York, 1965. 

8. Viskanta, R., “ Radiantion Transfer and Interaction of Convection with Radiantion 
Heat Transfer”, Advances in Heat Transfer , Vol. 3, Academic Press, New York, 1966. 

9. Tien, C.L., “Thermal Radiantion Properties of Gases”, Advances in Heat Transfer , 
Vol. 5, Academic Press, New York, 1968. 

10. Cess, R.D. and Tiwari, S.N., “Infrared Radiantive Energy Transfer in Gases”, 
Advances in Heat Transfer, Vol. 8, Academic Press, New York, 1972. 




70 


11. Edwards, D.K., “Molecular Gas Band Radiantion”, Advances in Heat Transfer , Vol. 
12, Academic Press, New York, 1976. 

12. Tiwari, S.N., “Band Models and Correlations for Infrared Radiantion”, Radian tive 
Transfer and Thermal Control (Progress in Astronautics and Aeronautics) , Vol. 49, 
American Institute of Aeronautics and Asrtonautics, New York, 1976. 

13. Tiwari, S.N., “Models for Infrared Atmospheric Radiantion”, Advances in Geo- 
physics , Vol. 20, Academic Press, New York, 1978. 

14. Tiwari, S.N., “Models for Infrared Atmospheric Radiantion”, Advances in Geo- 
physics , Vol. 20, Academic Press, New York, 1978. 

15. Plass, G.N., “Models for Spectral Band Absorption”, Journal of the Optical Society 
of America , Vol. 8, No. 10, October 1958, pp. 690-703. 

16. Plass, G.N., “Useful Representation for Measurments of Spectral Band Absorption”, 
Journal of the Optical Society of America , Vol. 50, No. 9, September 1960, pp. 
868-875. 

17. Wyatt, P.J., Stull, V.R., and Plass, G.N., “Quasi-Random Model of Band Absorp- 
tion”, Journal of the Optical Society of America , Vol. 52, No. 11, November 1962, 
pp. 1209-1217. 

18. Edwards, D.K. and Menard, W.A., “Comparison of Methods for Correlation of Total 
Band Absorption”, Applied Optics , Vol. 3, No. 5, May 1964, pp. 621-625. 

19. Edwards, D.K., Glassen, L.K., Hauser, W.C., and Tucher, J.S., “ Radiation Heat 
Transfer in Nonisothermal Nongray Gases”, Journal of Heat Transfer , Vol. 89C, No. 
3, 1967, pp. 219-229. 

20. Edwards, D.K. and Balakrishnan, A., “Slab Band Absorptance for Molecular Gas 


Radiantion”, Journal of the Quantitative Spectroscopy and Radiative Transfer , Vol. 
12, 1972, pp. 1379-1387. 



71 


21. Edwards, D.K. and Balakrishnan, A., “Thermal Radiantion by Combustion Gases”, 
International Journal of Heat and mass Transfer , Vol. 17, No. 1, January 1974, pp. 
155-158. 

22. Goody, R.M. and Belton, M.J.S., “ Radiantive Relaxation Times for Mars(A 
Discussion of Martian Atmospheric Dynamics)”, Planetary and Space Science , Vol. 
lr, No. 2, 1967, pp. 247-256. 

23. Tien, C.L. and Ling, G.R., “On a Simple Correlation for Total Band Absorptance 
of Radianting Gases”, International Journal of Heat and mass Transfer , Vol. 12, No. 
9, September 1969, pp. 1179-1181. 

24. Tiwari, S.N., “Applications of Infrared Band Model Correlations to Nongray 
Radiantion”, International Journal of Heat and mass Transfer , Vol. 20, No. 7, July 
1977, pp. 741-751. 

25. Ussikin, C. and Sparrow, E.M., “Thermal Radiantion Between Parallel Plates 
Separated by an Absorbing-Emitting Non-Isothermal Gas”, International Journal of 
Heat and mass Transfer , Vol. 1, 1960, pp. 28-36. 

26. Sparrow, E.M., Ussikin, C.M., and Hubbard, H.A., “ Radiantion Heat Transfer in 
a spherical Enclosure Containing a Participating, Heat-Generating Gas”, Journal of 
Heat Transfer , Vol. 82 C, 1961, pp. 197-206. 

27. Viskanta, R., “Interaction of Heat Transfer by Conduction, Convection, and 
Radiantion in Radianting Fluids”, Journal of Heat Transfer , Vol. 85 C, 1963, pp. 
318-328. 

28. Cess, R.D. and Tiwari, S.N., “Heat Transfer to Laminar Flow of an Absorbing- 
Emitting Gas between Parallel Plates”, Heat and Mass Transfer-USSR, Vol. 1, May 
1968, pp. 229-283. 



72 


29. Tiwari, S.N., “ Radiantive Interaction in Transient Energy Transfer in Gaseous Sys- 
tem , Dept, of Mechanical Engineering and Mechanics, Old Dominion University, 
Norfolk, Virginia, Progress Report NAG- 1-423, December 1985; also NASA-CR- 
176644 NAS 1.26:176644, December 1985. 

30. Tiwari, S.N. and Singh, D.J., “Interaction of Transient Radiantion in Fully Developed 
Laminar Flows”, AIAA Paper 87-1521, June 1987; also in Applied Scientific 
Research , December 1989. 

31. Einstein, T.H., Radiant Heat Transfer in Absorbing Gases enclosed in a Circular 
Pipe with Conduction, Gas Flow and Internal Heat Generation”, NASA TR R-156, 
1963. 

32. Nichols, L.D., “ Temperature Profile in the Entrance Region of an Annular Passage 
Considering the Effects of Turbulent Convection and Radiantion”, International 
Journal of Heat and mass Transfer , Vol. 8, 1965, pp. 589-607. 

33. deSoto, S. and Edwards, D.K., “ Radiantive Emission and Absorption in Nonisother- 
mal Nongray Gases in Tubes”, Proceedings 1965 Heat Transfer and Fluid Mechanics 
Institute , Stanford University Press, Stanford, Calif., 1965, pp. 358-372. 

34. deSoto, S., Coupled Radiantion, Conduction and Convection in Entrance Region 
Row”. International Journal of Heat and mass Transfer , Vol. 11, 1968, pp. 39-53. 

35. Pearce, B.E. and Emery, A.F., “Heat Transfer by Thermal Radiantion and Laminar 
Forced Convection to an Absorbing Fluid in the Entry Region of a Pipe”, Journal 
of Heat Transfer . Vol. 92, 1970, pp. 221-230 


APPENDICES 


73 



APPENDIX A 


DERIVATION OF CONSTANTS FOR CIRCULAR DUCT 

To determine constants in Eq.(5.6), Eq.(5.8) is evaluated at any two convinient 
locations. To avoid excessive writing, the following notations are used 


h, = bu 0 , 
c _ - - 1 


b\i@ 


(A. 1) 


where 




Rewriting governing Equation for circular duct using the quantities defined in (A.l) 
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Assuming solution of the form (A quartic solution in £) 


HO = a 1 (0-0 + a 2 (0-i) 

0'(O= 2Ai( + 4A 2 0 (A. 3) 
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By defining u = (£ - £') for the first and the fourth integrals, u = (£' - O 

for the second and the fifth integrals, and u = (£ 4 - £' — 2 £ sin 7 ) for the third and 

the sixth integrals in Eq.(A.4) and changing the limits respectively. And substituting for 
8 ^ j from (A. 3) into Eq.(A.2) and rearranging, equation in terms of changed variables is 
obtained. The integrals with the changed variables are solved individually and then added 
up as defined in Eqs. (5.9) and (5.10). The procedure adopted to evaluate constants ai 
and a 2 is similar to one applied in Ref. 39. 
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APPENDIX B 


ADDITIONAL RESULTS FOR LAMINAR FLOW OF NONGRAY 
GAS THROUGH PARALLEL PLATE PLATE GEOMETRY 

Extensive results obtained for laminar flow of nongray gas between two parallel 
black plates are presented here for reference purposes. The results are for species CO, 
OH, CO 2 . Results for temprature variations across the duct as well as variation of bulk 
temprature with plate spacing and pressure are presented. All the results show similar 
trend as discussed for Figs. 6.1-6.12. 
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Figure B.4 Variation of bulk temperature with plate spacing for CO 2 . 
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Figure B.5 Variation of bulk temperature with plate spacing for CO 2 ; 
P = 1 atm. 
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Figure B.8 Variation of 
P = 1 atm. 
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Figure B.l 1 Variation of bulk temperature with plate spacing for OH; 
P = 1 atm. 
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APPENDIX C 

ADDITIONAL RESULTS FOR LAMINAR FLOW 
OF NONGRAY GAS THROUGH CIRCULAR DUCT 


Extensive results obtained for laminar flow of nongray gas black circular duct are 
presented here for reference purposes. The results are for species CO, OH, C0 2 . Results 
for temprature variations across the duct as well as variation of bulk temprature with 
plate spacing and pressure are presented. All the results show similar trend as discussed 
for Figs. 6.13—6.22. 
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Figure C.5 Variation of bulk temperature with pressure for CO 2 ; 
T w = 1000 K and r 0 = 5 cm. 
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Figure C.6 Variation of bulk temperature with radius for CO. 
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Figure C.8 Variation of bulk temperature with press' 
T w = 1000 K and r 0 = 5 cm. 
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APPENDIX D 

PROGRAM LISTING FOR LAMINAR FLOW OF 
NONGRAY GAS THROUGH PARALLEL PLATE DUCT 

The listing of the FORTRAN program used for numerical procedure is given here. 
The program is coded on UNIX based SUN 386i. Input file consists of all the band 
parameters for different species[29]. The program uses romberg cautious integration 
routine for numerical integration purposes. The listing of romberg integration is not 
given here, but is available in IMSL libraries. 


program PLATE 


C kfb Thermal Conductivity of any Gas 

real l,kfb 

C n No. of Bands for the Gas 


parameter (n=2) 

external fpkO, fpkl, fpk2, fpk3, fpk4 
dimension fi (n) ,bi (n) ,hi (n) , eps (3) 
common /func/f,b 
open (5, file=' INPUT' ) 
rewind (5) 

open (6, file=' OUTPUT' ) 


rewind (6) 

C eps(l) Absolute Error Desired 

C eps (2) Relative Error Desired 

eps (1) = 1 . Oe-6 
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eps (2) = 1 . Oe-6 
do 333 i = 1, 656 

read (5, *) t,pr, 1, (fi (k) , k=l, n) , (bi (k) , k=l , n) , 

+ (hi (k) , k=l,n) , kfb 
am = 1/kfb 
hsum = 0.0 
suml = 0.0 
sum2 = 0.0 
sum3 = 0.0 
sum4 = 0.0 
do 111 k = 1, n 
f = fi(k) 
b = bi (k) 
h - hi (k) 
hsum= hsum+h 
cu = l./b 
cus = cu*cu 
cue = cu*cus 
cuq = cu*cuc 
x = 0.0 

• C Numerical Integration 

C Cautious Romberg Integration 

call cadre (x, b, fpkl, eps, 0, rrl, irl) 
call cadre (x, b, fpk2, eps, 0, rr2, ir2) 
call cadre (x, b, fpk3, eps, 0, rr3, ir3) 
call cadre (x,b, fpk4, eps, 0, rr4, ir4) 
suml = suml+ (h* ( (cu*rrl) - (2 . *cuc*rr3) + (cuq*rr4) ) ) 
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sum2 = sum2+ (h* ( (cus*rr2) - (2 . *cuc*rr3) + (cuq*rr4) ) ) 
p = b/4 . 
q = (3./4.)*b 

call cadre (p, q, fpkO, eps, 0, ssl, isl) 
call cadre (x,p, fpkl, eps,0,ss2,is2) 
call cadre (x, q, fpkl, eps, 0,ss3, is3) 
call cadre (p, q, fpk2, eps, 0, ss4, is4) 
call cadre (x,p,fpk3,eps,0,ss5,is5) 
call cadre (p, q, fpk3, eps, 0,ss6, is6) 
call cadre (p, q, fpk4, eps, 0, ss7, is7) 
brkl = (S7./256. ) *ssl 
brk2 = (11 . /16 . ) *cu* (ss2 + ss3) 
brk3 = ( 9 . /8 . ) *cus*ss4 
brk4 = cue* (ss5+ss6) 
brk5 = cuq*ss7 

sum3 = sum3+ (h* (brkl+brk2-brk3-brk4+brk5) ) 

brsl = (9./256.) *ssl 

brs2 = (3./16. ) *cu* (ss2 + ss3) 

brs3 = (1 . /8 . ) *cus*ss4 

brs5 = cuq*ss7 

sum4 = sum4+ (h* (brsl+brs2-brs3-brs4+brs5) ) 

111 continue 

alphal = 1 . 0+ (am*suml ) 
alpha2 = am*sum2 
alpha3 = (11 . / 1 6 . ) + (am*sum3) 
alpha4 = ( 3 . /1 6 . ) + (am*sum4 ) 

constt = 16 . * ( (alphal*alpha4) - (alpha2*alpha3) ) 
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constl = (11 . *alpha2) - (16 . *alpha4) 
const2 = (16 . * alpha 3) - (11 . * alpha 1) 
al = const 1 /constt 
a2 = const2 /constt 
bulk = -((17.*al)+(3.*a2))/70. 
call lpl (am,hsum, coni, con2,bulkl) 
write (6, *) t, pr, 1, bulk, bulkl 
else 

continue 

endif 

333 continue 
end 

function fpkO (u) 
common /func/f,b 

den = (f* ( (u*u) + (2.*u)+2.)+u) * (u+(2.*f) ) 

aud = (f* ( (u*u) + (4.*u*f) + (4.*f) ) ) /den 

fpkO= aud 

return 

end 

function fpkl (u) 
common /func/f,b 

den = (f * ( (u*u) + (2 . *u) +2 . ) +u) * (u+ (2 . *f ) ) 

aud = (f * ( (u*u) + (4 . *u*f) + (4 . *f) ) ) /den 

fpkl= u*aud 

return 

end 

function fpk2 (u) 
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common /func/f,b 

den = (f* ( (u*u) + ( 2 .*u) +2 . ) +u) * (u+ (2 . *f ) ) 

aud = (f * ( (u*u) + (4 . *u*f ) + (4 . *f ) ) ) /den 

fpk2= u*u*aud 

return 

end 

function fpk3 (u) 
common /func/f,b 

den = (f * ( (u*u) + (2 . *u) +2 . ) +u) * (u+ (2 . *f ) ) 

aud = (f * ( (u*u) + (4 . *u*f ) + (4 . *f ) ) ) /den 

fpk3= u*u*u*aud 

return 

end 

function fpk4 (u) 
common /func/f,b 

den = (f* ( (u*u) + (2 . *u) +2.) +u) * (u+(2.*f) ) 

aud = (f * ( (u*u) + (4 . *u*f ) + (4 . *f ) ) ) /den 

fpk4= u*u*u*u*aud 

return 

end 

subroutine lpl ( am, bt , coni , con2 , bulkl ) 
em = am*ht 

etal = l.+((7./12.)*era) 

eta2 = (1./12 .) *em 

eta3 = (ll./16.)+(0.583154559*em) 

eta4 = (3 . /16 . ) + (0 . 127l64755*em) 

phfc = 16 . * ( (etal*eta4) — (eta2*eta3) ) 
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phil = (11 . *eta2) - (16 . *eta4) 
phi2 = (16. *eta3) - (11 . *etal) 
coni = phil/phic 
con2 = phi2/phic 

bulkl= -((17.*conl)+(3. *con2) ) /70 . 

return 

end 


106 


APPENDIX E 

PROGRAM LISTING FOR LAMINAR FLOW OF 
NONGRAY GAS THROUGH CIRCULAR DUCT 
The listing of the FORTRAN program used for numerical procedure is given here. 
The program is coded on UNIX based SUN 386i. Input file consists of all the band param- 
eters for different species[29]. The program uses Multidimensional quadrature integration 
routine using Gauss-Kronrod rule (TWODQ), for numerical integration purposes. The 
subroutine is available in IMSL library (Chap. 4, Vol. 2.). 

program TUBE 
real kfb 

C kfb Thermal Conductivity of the Gas 

parameter (n=l,pi=3. 1415926535898) 

C n No. of Bands in any Gas 

external frll, frl2, fr21, fr22, fr31, fr32 
external o,q, r,s 

dimension f i (n) , bi (n) , hi (n) , zi (2) 
dimension zz (2) , gamal (2) , gama2 (2) 
common /func/f,b,z 

C INPUT Input Data File with All Properties for 

C the Gas under consideration. 

open <5, file*' INPUT' ) 
rewind (5) 

open (6, file*' OUTPUT' ) 
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rewind (6) 

C ea Absolute Error Desired 

C er Relative Error Desired 

ea = 1.0e-4 
er = 1 . Oe-4 

C x Lower Limit of Integration for Outer Integral 

C y Upper Limit of Integration for Outer Integral 

x = 0.0 
y = (pi/2.0) 

C zi(l) & zi(2) Two Values of xi . 

zi (1) = 0.5 
zi (2) = 1.0 
do 333 i = l,imax 

read (5, *) t,p, 1, (fi (k) , k=l,n) , (bi (k) , k=l,n) , 

+ (hi (k) ,k=l,n) ,kfb 
am = (3.2*1) / (pi*kfb) 
do 222 k = 1,2 
z = zi (k) 
xz = 2 . *z 
yz = 4 . * (z**3 . ) 
zz (k) = (z**3 . ) - (2 . *z) 
brkl =0.0 
brk2 =0.0 
hsum = 0.0 
do 111 j = l,n 
f = fi ( j) 
b = bi ( j) 
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h = hi ( j) 
hsum= hsum+h 


C Numerical Integration 

C Multidimensional Quadrature Technique using 

C Gauss-Kronrod rule. 


call twodq ( f r 1 1 , x, y , o, q, ea, er, 1 , rrl 1 , erl 1 ) 
call twodq (frl 2 ,x,y,o,q,ea,er,l,rrl 2 ,erl 2 ) 
call twodq (fr21, x, y , o, r,ea,er, 1, rr2 1 , er21 ) 
call twodq (fr22, x, y, o, r, ea, er, 1, rr22, er22) 
call twodq (fr 31 , x, y, q, s,ea,er, I,rr31,er31) 
call twodq (fr32,x,y,q,s,ea,er,l,rr32,er32) 
brkl = brkl+ (h* (rrll-rr21+rr31) ) 
brk2 = brk2+ (h* (rrl2-rr22+rr32) ) 

111 continue 

gamal (k) = (am*brkl)-xz 

gama2(k) = (am*brk2) -yz 
222 continue 

alphl = ( (zz ( 1 ) *gama2 (2) ) - (zz (2) *gama2 (1) ) ) 

alph2 = ( (zz ( 2 ) *gamal ( 1 ) ) - (zz (1) *gamal (2) ) ) 

alphc = ( (gamal ( 1 ) *gama 2 ( 2 ) )-(gama2 (1) *gamal (2) ) ) 

al = alphl/alphc 
a2 = alph2/alphc 
em = hsum*am 

call lpl (em, coni, con2, bulkl) 

bulk = - ( ( (-2 . /3 . ) *al) + ( (-5./6.)*a2) ) 

bulkl = - ( ( (- 2 . / 3 . ) *conl) + ( (-5 . /6 . ) *con2) ) 
write(6, *)t,p,r, bulk, bulkl 
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if (r . eq. 5 . 0) then 
do 11 ii = 1, 11 
xi = 0 . l*f loat (ii-1) 

theta = -(al*((xi**2.)-l.)+a2*((xi**4.)-l.)) 
thetal= -(coni* ( (xi**2.)-l.) +con2* ( (xi**4 . ) -1 . ) ) 
write ( 6, *) t , p, xi, theta, thetal 
11 continue 
else 

continue 

endif 

333 continue 
end 

function frll(x,u) 
common /func/f,b,z 

den = (f*((u*u) + (2.*u)+2.)+u)*(u+(2.*f)) 
aud = (f * ( (u*u) + (4 . *u*f ) + (4 . *f ) ) ) /den 
frll= ( (z- (cos (x) *u/b) ) **2.-1.) *cos (x) *aud 
return 
end 

function frl2(x,u) 
common /func/f,b,z 

den = (f* ( (u*u) + (2.*u) +2.)+u) * (u+ (2.*f) ) 
aud = (f * ( (u*u) + (4 . *u*f ) + (4 . *f ) ) ) /den 
frl2= ( (z-(cos (x) *u/b) ) **4.-1.) *cos (x) *aud 
return 
end 

function fr21(x,u) 


110 



common /func/f,b,z 

den = (f * ( (u*u) + (2 . *u) +2 . ) +u) * (u+ (2 . *f ) ) 
aud = (f* ( (u*u) + (4.*u*f) + (4.*f) ) ) /den 

f r21= ( (z+ (cos (x) *u/b) ) **2.-1 . ) *cos (x) *aud 
return 
end 

function fr22(x,u) 
common /func/f,b,z 

den = (f * ( (u*u) + ( 2 . *u) +2 . ) +u) * (u+ (2 . *f ) ) 

aud = (f * ( (u*u) + (4 . *u*f ) + (4 . *f ) ) ) /den 
fr22= ( (z+ (cos (x) *u/b) ) ** 4 . -1 . ) *cos (x) *aud 
return 
end 

function fr31(x,u) 
common /func/f,b,z 

den = (f*((u*u) + (2.*u)+2.)+u)*(u+(2.*f)) 

aud = (f*((u*u) + (4.*u*f) + (4.*f)))/den 

f r31= < (z* (2 . *sin (x) -1 . ) + (cos (x) *u/b) >**2.-1.) *cos (x) *aud 

return 

end 

function fr32(x,u) 
common /func/f,b,z 

den = (f * ( (u*u) + (2 . *u) +2 . ) +u) * (u+ (2 . *f ) ) 
aud = (f*((u*u) + (4.*u*f) + (4.*f)))/den 

f r32= { (z*(2.*sin (x) -1.) + (cos * u / b > ) **4 .-1 . ) *cos (x) *aud 

return 

end 
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function o{x) 

common /func/f,b,z 

o = 0.0 

return 

end 

function q(x) 

common /func/f,b,z 

q = b*z* ( 1 . -sin (x) ) /cos (x) 

return 

end 

function r (x) 
common /func/f,b,z 
r = b* ( 1 . -z) /cos (x) 
return 
end 

function s (x) 
common /func/f,b,z 

s = b* (1 . +z- (2 . *z*sin (x) ) ) /cos (x) 

return 

end 
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